
# Replication Files

# "Foreign Policy Failures and the Attractiveness of Great Powers"
# by Rachel Myrick and William Marble, BJPolS 

# Main analysis script of Gallup World Poll data


library(tidyverse)
library(fixest)
library(lubridate)
library(rio)
library(tinytable)
library(rdrobust)
library(modelsummary)
# hypotheses test for seemingly unrelated regressions
if (!require(vcovSUR)) {
  devtools::install_github("kylebutts/vcovSUR")
  require(vcovSUR)
}
library(ggsignif)
library(texreg)
options(modelsummary_format_numeric_latex = "plain")



theme_set(theme_classic() + 
            theme(axis.text = element_text(size = 14), 
                  title = element_text(size = 16),
                  strip.text = element_text(size = 16),
                  legend.text = element_text(size = 14)))

treatment_date <- ymd("2021-08-15")


# Read GWP data -----------------------------------------------------------

# Due to licensing agreements, we cannot distribute this raw dataset.
# See https://guides.library.upenn.edu/gallup
# https://www.gallup.com/analytics/318875/global-research.aspx

# For data description
origdat <- import("Data/Gallup_World_Poll_Wave_11_17_021723/Gallup_World_Poll_Wave_11_17_021723.sav")

dat <- origdat %>%
  mutate(
    date_interview = WP4,
    country = WP5,
    caseid = WPID,
    gender = WP1219,
    hhsize = WP12,
    gender2 = WP1219,
    age = WP1220,
    marital = WP1223,
    educ2 = WP3117,
    borncountry = WP4657,
    leader_usa = WP151,
    leader_eu = WP1423,
    leader_germany = WP153,
    leader_france = WP154,
    leader_russia = WP155,
    leader_china = WP156,
    leader_japan = WP157,
    leader_india = WP162,
    leader_turkey = WP167,
    leader_iran = WP187,
  ) 



## Covariates to test for balance ------------------------------------------
dat <- dat %>% 
  mutate(age = as.numeric(age),
         male = as.integer(gender == 1),
         educ_secondary = case_when(
           educ2 %in% c(2,3) ~ 1L,
           educ2 == 1 ~ 0L,
           educ2 >= 4 ~ NA
           ),
         educ = case_when(
           educ2 == 1 ~ "Primary school",
           educ2 == 2 ~ "Secondary school",
           educ2 == 3 ~ "Postsecondary school",
           TRUE ~ NA
         ),
         inc_quint_1 = as.numeric(INCOME_5 == 1),
         inc_quint_2 = as.numeric(INCOME_5 == 2),
         inc_quint_3 = as.numeric(INCOME_5 == 3),
         inc_quint_4 = as.numeric(INCOME_5 == 4),
         inc_quint_5 = as.numeric(INCOME_5 == 5),
         
         hhsize = ifelse(hhsize >= 30, 30, hhsize) # top code hhsize at 30... there are a few of 98+ but seems wrong
         ) %>% 
  mutate(educ = factor(educ, c("Primary school", "Secondary school", "Postsecondary school")))
covs <- c("age", "male", "educ", 
          "inc_quint_1", "inc_quint_2",  "inc_quint_3", "inc_quint_4", "inc_quint_5",
          "hhsize")


## Main variables ------------------------------------------------------------
dat <- dat %>% 
  mutate(across(starts_with("leader"), ~ {
    case_when(
      .x == 1 ~ "Approve",
      .x == 2 ~ "Disapprove",
      TRUE ~ NA_character_
    )
  }),
  after_withdrawal = as.integer(date_interview >= treatment_date),
  days_withdrawal  = date_interview - treatment_date
  )

dat <- dat %>% 
  mutate(COUNTRYNEW = ifelse(COUNTRYNEW == "Bosnia Herzegovina", "Bosnia and Herzegovina", COUNTRYNEW))

# create binary indicators for approval
dat <- dat %>% 
  mutate(across(starts_with("leader"), list(bin = ~ as.numeric(.x == "Approve"))))

## all data (w cleaned names)
alldat <- dat

## Subset to responses during Biden administration that had an answer 
# on the US leadership question
dat <- dat %>% 
  filter(!is.na(leader_usa),
         date_interview >= ymd("2021-03-01"))







# Plot full time series ---------------------------------------------------


## Plot approval before/after
ggplot(dat) + 
  aes(x = date_interview, 
      y = leader_usa_bin, 
      weight = HHWEIGHT2, 
      group = factor(after_withdrawal)) + 
  geom_vline(xintercept = treatment_date) + 
  geom_smooth() + 
  scale_x_date(date_breaks = "6 months", 
               date_labels = "%b %Y") + 
  scale_y_continuous(labels = scales::label_percent()) + 
  labs(x = NULL, y = "Approval of U.S. Leadership")
ggsave("Figures/US-approval-time-series.pdf", width = 6, height = 4)


## Residualize by country
mod <- lm(leader_usa_bin ~ COUNTRYNEW,
          weights = HHWEIGHT2,
          data = dat)
dat$leader_usa_resid <-  mod$residuals

ggplot(dat %>% 
         filter(date_interview <= ymd("2023-08-15"))) + 
  aes(x = date_interview, 
      y = leader_usa_resid, 
      weight = HHWEIGHT2, 
      group = factor(after_withdrawal)) + 
  geom_vline(xintercept = treatment_date) + 
  stat_summary_bin(binwidth = 30) +
  geom_smooth(formula = y ~ poly(x, 4)) +
  annotate(x = treatment_date+10, 
           y = .1, 
           geom = "text", 
           label = "Afghanistan\nWithdrawal",
           hjust = 0) + 
  scale_x_date(date_breaks = "6 months", 
                expand = expansion(), 
               date_labels = "%b '%y") + 
  scale_y_continuous(labels = scales::label_percent()) + 
  labs(x = NULL, y = "Approval of U.S. Leadership\n(residualized by country)")
ggsave("Figures/US-approval-resid-time-series.pdf", width = 6, height = 4)






# Plot w/i bandwidth ------------------------------------------------------


# Subset to responses 30 days before/after Kabul fell
bw <- 30
datbw <- dat %>% 
  filter(date_interview >= treatment_date - bw,
         date_interview <  treatment_date + bw) 

# Get codebook for +/- 30 days
labs <- lapply(names(datbw), function(x) {
  if (x != toupper(x)) {
    return(NULL)
  } else if (all(is.na(datbw[[x]]))) {
    NULL
  } else if (is.null(attr(datbw[[x]], "label"))) {
    return(NULL)
  } else {
    data.frame(var = x, label = attr(datbw[[x]], "label"))
  }
}) %>% 
  bind_rows()




## Get Sample Sizes by Country X Before/After ----------



## Keep countries w/ at least 25 interviews both before/after
sampsizes <- datbw %>% 
  filter(!is.na(leader_usa)) %>% 
  count(COUNTRYNEW, factor(after_withdrawal), .drop = FALSE) %>% 
  group_by(COUNTRYNEW) %>% 
  mutate(min = min(n)) %>% 
  filter(min >= 25)
tokeep <- unique(sampsizes$COUNTRYNEW)
tokeep <- tokeep[-which(tokeep == "Afghanistan")]
datbw <- datbw %>% 
  filter(COUNTRYNEW %in% tokeep)
sampsizes <- sampsizes %>% 
  filter(COUNTRYNEW %in% datbw$COUNTRYNEW)

# full sample size
ntot <- prettyNum(sum(sampsizes$n), big.mark = ",")
cat(ntot, file = "Tables/n-total.tex")

## Make Table showing sample sizes
sstab <- sampsizes %>% 
  select(-min) %>% 
  pivot_wider(values_from = n, names_from = `factor(after_withdrawal)`) %>% 
  rename(Country = COUNTRYNEW, `Pre-Withdrawal` = `0`, `Post-Withdrawal` = `1`) %>% 
  mutate(Total = `Pre-Withdrawal`+ `Post-Withdrawal`) %>% 
  arrange(desc(Total), Country)
n_countries <- length(unique(sstab$Country))
cat(n_countries, file = "Tables/n-countries.tex")
cat(sort(unique(sstab$Country)), sep = ", ", file = "Tables/country-list.tex")
add <- sstab %>% 
  ungroup() %>% 
  select(-Country) %>% 
  summarise(across(everything(), sum)) %>% 
  mutate(Country = "Total")
sstab <- bind_rows(sstab, add)

sstab %>% 
  # mutate(across(2:4, ~ prettyNum(.x, big.mark = ","))) %>% 
  xtable::xtable() %>% 
  print(include.rownames = FALSE, 
        format.args =  list(big.mark = ","),
        sanitize.colnames.function = function(x) {paste('{\\textbf{',x,'}}', sep ='')},
        booktabs = TRUE, 
        floating = FALSE,
        add.to.row = list(pos = list(-1, nrow(sstab)), 
                          command = c("\\toprule", "\\bottomrule")),
        hline.after = c(0, nrow(sstab)-1)) %>% 
  cat(file = "Tables/Sample-Sizes.tex")



## Within-country averages ---------------------------------

# Calculate approval of US before/after
outlist <- list()
for (cntry in unique(datbw$COUNTRYNEW)){
  tmp <- filter(datbw, COUNTRYNEW == cntry)
  before <- filter(tmp, after_withdrawal == 0)
  after <- filter(tmp, after_withdrawal == 1)
  
  test <- weights::wtd.t.test(
    x = (before$leader_usa_bin),
    weight = before$HHWEIGHT2,
    y = (after$leader_usa_bin),
    weighty = after$HHWEIGHT2,
    samedata = FALSE
  )
  
  out <- tibble(country = cntry,
                approve_before = test$additional["Mean.x"],
                approve_after = test$additional["Mean.y"],
                se = test$additional["Std. Err"],
                p = test$coefficients["p.value"],
                n = nrow(before) + nrow(after)) 
  
  outlist[[cntry]] <- out
}
xtab <- bind_rows(outlist)



ggplot(xtab %>% 
         mutate(country = case_when(
           p >= .1 ~ country,
           p < .1 & p >= .05 ~ paste0( "*", country),
           p < .05 & p >= .01 ~ paste0("**", country),
           p < .01 ~ paste0("***", country)
         )) %>% 
         mutate(country = paste0(country, " (", prettyNum(n, big.mark = ","), ")"))) + 
  aes(x = approve_before,
      xend = approve_after,
      y = fct_reorder(country, approve_after),
      yend = fct_reorder(country, approve_after)) + 
  geom_segment(arrow = arrow(length = unit(.1, "inches"))) + 
  labs(x = "Approval of U.S. Leadership in 30 Days\nBefore/After August 15, 2021",
       y = NULL) +
  coord_cartesian(clip = "off",
                  ylim = c(0, 26)) + 
  scale_x_continuous(breaks = seq(0, 1, .1), labels = scales::label_percent(),
                     limits = c(.1,1)) +  
  theme(panel.grid.major.y = element_line()) + 
  annotation_custom(
    grob = grid::textGrob(label = "Country (Sample Size)", hjust = .63, gp = grid::gpar(cex = 1.25)),
    ymin = 25.75,
    ymax = 25.75,
    xmin = -0.1,
    xmax = -0.1)
ggsave("Figures/US-approval-change-by-country.pdf",
       width = 8, height = 6)


## Overall average in estimation sample ---------------------------

pct <- 100*weighted.mean(datbw$leader_usa_bin, datbw$HHWEIGHT2)
cat(sprintf("%s\\%%", round(pct, 1)), file = "Tables/overall-mean.tex")



# Regression --------------------------------------------------------------




base <- feols(leader_usa_bin ~ after_withdrawal,
            cluster = ~ COUNTRYNEW,
            weights = ~ HHWEIGHT2,
            data = datbw)

basefe <- feols(leader_usa_bin ~ after_withdrawal | COUNTRYNEW,
                cluster = ~ COUNTRYNEW,
                weights = ~ HHWEIGHT2,
                data = datbw)

basefecov <- feols(leader_usa_bin ~ after_withdrawal | COUNTRYNEW[age] + 
                     COUNTRYNEW^educ + COUNTRYNEW^gender2,
                  cluster = ~ COUNTRYNEW,
                  weights = ~ HHWEIGHT2,
                  data = datbw)
rd <- feols(leader_usa_bin ~ days_withdrawal * after_withdrawal, 
            cluster = ~ COUNTRYNEW,
            weights = ~HHWEIGHT2,
            data = datbw)
rdfe <- feols(leader_usa_bin ~ days_withdrawal * after_withdrawal | COUNTRYNEW, 
              cluster = ~ COUNTRYNEW,
              weights = ~HHWEIGHT2,
              data = datbw)
rdfecov <- feols(leader_usa_bin ~ days_withdrawal * after_withdrawal | COUNTRYNEW[age] + 
                   COUNTRYNEW^educ + COUNTRYNEW^gender2, 
                  cluster = ~ COUNTRYNEW,
                  weights = ~HHWEIGHT2,
                  data = datbw)

glance_custom.fixest <- function(x, ...) {
  # tibble::tibble(`N_Countries` = paste(x$fixef_sizes[COUNTRYNEW], collapse = " + "))
}




modelsummary(list(base, basefe, basefecov, rd, rdfe, rdfecov), 
             stars = TRUE)

modelsummary(list(base, basefe, basefecov, rd, rdfe, rdfecov), 
             stars = c("$^*$" = .1, "$^{**}$" = .05, "$^{***}$" = .01), 
             coef_map = c("after_withdrawal" = "After Withdrawal",
                          "days_withdrawal" = "Time",
                          "days_withdrawal:after_withdrawal" = "Time $\\times$ After Withdrawal"),
             booktabs = TRUE, 
             gof_omit = c("Std.Errors|AIC|BIC|R2 Adj."),
             output = "latex_tabular",
             escape = FALSE,
             gof_map = tribble(~ raw, ~ clean, ~ fmt, 
                               "nobs", "$N$", function(.x) prettyNum(.x, big.mark = ","),
                               "FE: COUNTRYNEW", "Country FE", modelsummary::fmt_statistic),
             add_rows = tribble(~ term, ~ base, ~ basefe, ~ basefecov, ~ rd, ~ rdfe, ~ rdfecov,
                                "Ind. Covariates $\\times$ Country FE", "", "", "X", "", "", "X")
             ) %>% 
  format_tt(replace = list("$\\checkmark$" = "X")) %>% 
  save_tt("Tables/RD-results.tex", overwrite = TRUE)
  




# Leave one country out at a time -----------------------------------------

res <- list()
countries <- unique(datbw$COUNTRYNEW)
for (cntry in countries) {
  subs <- datbw %>% filter(COUNTRYNEW != cntry)
  rdfecov <- feols(leader_usa_bin ~ days_withdrawal * after_withdrawal | COUNTRYNEW[age] + 
                     COUNTRYNEW^educ + COUNTRYNEW^gender2, 
                   cluster = ~ COUNTRYNEW,
                   weights = ~HHWEIGHT2,
                   data = subs)
  res[[cntry]] <- tidy(rdfecov) %>% mutate(country_left_out = cntry)
}
res <- bind_rows(res)
ggplot(res %>% filter(term == "after_withdrawal")) + 
  aes(y = fct_rev(country_left_out),
      x = estimate, 
      xmin = estimate - 1.96 * std.error, xmax = estimate + 1.96 * std.error) +
  geom_vline(xintercept = 0, lty = 2) + 
  geom_point() + 
  geom_errorbarh(height = 0) + 
  guides(colour = "none") + 
  scale_x_continuous(breaks = seq(-.1, .1, .02), limits = c(-0.1, 0.05)) +
  labs(x = "RD Estimate", y = "Country Left Out of Estimation")
ggsave("Figures/RD-leave-out-countries.pdf", width=8, height=6)







# Varying BW --------------------------------------------------------------

bws <- 5:100
results <- list()
i = 1
for (bw in bws){
  tmpdat <- dat %>% 
    filter(!is.na(leader_usa), 
           date_interview >= treatment_date - bw,
           date_interview <  treatment_date + bw)
  
  
  
  ## Keep countries w/ at least 25 interviews both before/after
  sampsizes <- tmpdat %>% 
    filter(!is.na(leader_usa)) %>% 
    count(COUNTRYNEW, factor(after_withdrawal), .drop = FALSE) %>% 
    group_by(COUNTRYNEW) %>% 
    mutate(min = min(n)) %>% 
    filter(min >= 25)
  tokeep <- unique(sampsizes$COUNTRYNEW)
  tokeep <- tokeep[-which(tokeep == "Afghanistan")]
  tmpdat <- tmpdat %>% 
    filter(COUNTRYNEW %in% tokeep)
  sampsizes <- sampsizes %>% 
    filter(COUNTRYNEW %in% tmpdat$COUNTRYNEW)
  # end(warning)
  
  
  ## Fit models
  mod1 <- feols(leader_usa_bin ~ days_withdrawal * after_withdrawal | COUNTRYNEW, 
                cluster = ~ COUNTRYNEW,
                weights = ~HHWEIGHT2,
                data = tmpdat)
  nc1 <- length(attr(mod1$fixef_id$COUNTRYNEW, "fixef_names"))
  
  
  # use i() syntax for age covariate w/ varying slopes rather than COUNTRYNEW[age] 
  # due to convergence issues with the latter.
  mod2 <- feols(leader_usa_bin ~ days_withdrawal * after_withdrawal + i(COUNTRYNEW, age) |
                  COUNTRYNEW + COUNTRYNEW^educ + COUNTRYNEW^gender2, 
               cluster = ~ COUNTRYNEW,
               weights = ~HHWEIGHT2,
               data = tmpdat)
  nc2 <- length(attr(mod2$fixef_id$COUNTRYNEW, "fixef_names"))
  
  
  out <- bind_rows(tidy(mod1) %>% mutate(model = "basic") %>% mutate(n_countries = nc1),
                   tidy(mod2) %>% mutate(model = "controls") %>% mutate(n_countries = nc2)) %>% 
    filter(term == "after_withdrawal") %>% 
    mutate(bw = bw) 
  
  results[[i]] <- out
  i = i+1
}
out <- bind_rows(results)

# include number of countries in the plot
ncs <- out %>% 
  distinct(bw, n_countries) %>% 
  filter(bw %% 10 == 0)

ggplot(out) + 
  aes(x = bw, 
      y = estimate, 
      ymin = estimate - 1.96 * std.error, 
      ymax = estimate + 1.96 * std.error,
      colour = model, 
      fill = model) + 
  geom_hline(yintercept = 0, lty = 3) +
  geom_ribbon(alpha = .05) + 
  geom_point() + 
  geom_text(data = ncs, mapping = aes(x = bw, label = n_countries), 
            y = 0.01, inherit.aes = FALSE) + 
  annotate(geom = "text", x = 7.5, y = .02, label = "N Countries:", hjust = 0) + 
  scale_x_continuous(breaks = seq(0, 100, 10)) + 
  scale_fill_manual(name = NULL, labels = c("w/ country FE", "w/ FEs + controls"), values = c("black", "blue")) + 
  scale_colour_manual(name = NULL, labels = c("w/ country FE", "w/ FEs + controls"), values = c("black", "blue")) + 
  labs(x = "Bandwidth (Days)", y = "RD Estimate") +
  theme(legend.position = c(.75, .15))
ggsave("Figures/RD-bandwith-plot.pdf", width = 6, height = 4)



# write smallest and largest effect sizes
cat(sprintf("%s", round(min(out$estimate[out$model == "controls"]), 3)), file = "Tables/min-effect.tex")
cat(sprintf("%s", round(max(out$estimate[out$model == "controls"]), 3)), file = "Tables/max-effect.tex")



# "Donut" regressions -------------------------------------------------------

## Exclude +/- 3 days from August 15

donutdat <- datbw %>% 
  filter(abs(days_withdrawal) > 2)

base <- feols(leader_usa_bin ~ after_withdrawal,
              cluster = ~ COUNTRYNEW,
              weights = ~ HHWEIGHT2,
              data = donutdat)
basefe <- feols(leader_usa_bin ~ after_withdrawal | COUNTRYNEW,
                cluster = ~ COUNTRYNEW,
                weights = ~ HHWEIGHT2,
                data = donutdat)
basefecov <- feols(leader_usa_bin ~ after_withdrawal | COUNTRYNEW[age] + 
                     COUNTRYNEW^educ2 + COUNTRYNEW^gender2,
                   cluster = ~ COUNTRYNEW,
                   weights = ~ HHWEIGHT2,
                   data = donutdat)
rd <- feols(leader_usa_bin ~ days_withdrawal * after_withdrawal, 
            cluster = ~ COUNTRYNEW,
            weights = ~HHWEIGHT2,
            data = donutdat)
rdfe <- feols(leader_usa_bin ~ days_withdrawal * after_withdrawal | COUNTRYNEW, 
              cluster = ~ COUNTRYNEW,
              weights = ~HHWEIGHT2,
              data = donutdat)
rdfecov <- feols(leader_usa_bin ~ days_withdrawal * after_withdrawal | COUNTRYNEW[age] + 
                   COUNTRYNEW^educ2 + COUNTRYNEW^gender2, 
                 cluster = ~ COUNTRYNEW,
                 weights = ~HHWEIGHT2,
                 data = donutdat)



modelsummary(list(base, basefe, basefecov, rd, rdfe, rdfecov), 
             stars = c("$^*$" = .1, "$^{**}$" = .05, "$^{***}$" = .01), 
             coef_map = c("after_withdrawal" = "After Withdrawal",
                          "days_withdrawal" = "Time",
                          "days_withdrawal:after_withdrawal" = "Time $\\times$ After Withdrawal"),
             booktabs = TRUE, 
             gof_omit = c("Std.Errors|AIC|BIC|R2 Adj."),
             output = "latex_tabular",
             escape = FALSE,
             gof_map = tribble(~ raw, ~ clean, ~ fmt, 
                               "nobs", "$N$", function(.x) prettyNum(.x, big.mark = ","),
                               "FE: COUNTRYNEW", "Country FE", modelsummary::fmt_statistic),
             add_rows = tribble(~ term, ~ base, ~ basefe, ~ basefecov, ~ rd, ~ rdfe, ~ rdfecov,
                                "Ind. Covariates $\\times$ Country FE", "", "", "X", "", "", "X")
) %>% 
  format_tt(replace = list("$\\checkmark$" = "X")) %>% 
  save_tt("Tables/RD-results-donut.tex", overwrite = TRUE)






# Zero-sum Reputation Analyses ----------------------------------------

## Raw approval before/after ------------------------------------------

# Calculate approval of Russia and China before/after US withdrawal
outlist_rc <- list()
for (cntry in unique(datbw$COUNTRYNEW)){
  tmp <- filter(datbw, COUNTRYNEW == cntry)
  before <- filter(tmp, after_withdrawal == 0)
  after <- filter(tmp, after_withdrawal == 1)
  
  test_russia <- weights::wtd.t.test(
    x = (before$leader_russia_bin),
    weight = before$HHWEIGHT2,
    y = (after$leader_russia_bin),
    weighty = after$HHWEIGHT2,
    samedata = FALSE
  )
  test_china <- weights::wtd.t.test(
    x = (before$leader_china_bin),
    weight = before$HHWEIGHT2,
    y = (after$leader_china_bin),
    weighty = after$HHWEIGHT2,
    samedata = FALSE
  )
  
  out_russia <- tibble(country = cntry,
                       approve_before = test_russia$additional["Mean.x"],
                       approve_after = test_russia$additional["Mean.y"],
                       se = test_russia$additional["Std. Err"],
                       p = test_russia$coefficients["p.value"],
                       n = nrow(before) + nrow(after), 
                       target = "Russia") 
  
  out_china <- tibble(country = cntry,
                      approve_before = test_china$additional["Mean.x"],
                      approve_after = test_china$additional["Mean.y"],
                      se = test_china$additional["Std. Err"],
                      p = test_china$coefficients["p.value"],
                      n = nrow(before) + nrow(after), 
                      target = "China") 
  
  outlist_rc[[cntry]] <- bind_rows(out_russia, out_china)
}
xtab_rc <- bind_rows(outlist_rc)



ggplot(xtab_rc %>% 
         mutate(country = paste0(country, " (", prettyNum(n, big.mark = ","), ")"))) +
  aes(x = approve_before,
      xend = approve_after,
      y = fct_reorder(country, approve_after),
      yend = fct_reorder(country, approve_after),
      colour = target) + 
  geom_segment(arrow = arrow(length = unit(.1, "inches")), lwd = 1.5, alpha = .7) + 
  labs(x = "Approval of Country Leadership in 30 Days\nBefore/After August 15, 2021",
       y = NULL) +
  coord_cartesian(clip = "off",
                  ylim = c(0, 26)) +
  scale_colour_manual(name = "Approval of...", values = c("firebrick", "black")) + 
  scale_x_continuous(breaks = seq(0, 1, .2), labels = scales::label_percent()) +  
  theme(panel.grid.major.y = element_line())  +
  annotation_custom(
    grob = grid::textGrob(label = "Country (Sample Size)", hjust = .63, gp = grid::gpar(cex = 1.25)),
    ymin = 25.75,
    ymax = 25.75,
    xmin = -0.2,
    xmax = -0.2)
ggsave("Figures/Russia-China-approval-change-by-country.pdf",
       width = 8, height = 6)


## Main RD Specifications ---------------------------------------------

vars <- c(
  "leader_usa_bin",
  "leader_eu_bin",
  "leader_germany_bin",
  "leader_france_bin",
  "leader_russia_bin",
  "leader_china_bin",
  "leader_japan_bin",
  "leader_india_bin",
  "leader_turkey_bin",
  "leader_iran_bin"
)

res <- list()
for (v in vars) {
  country <- tools::toTitleCase(gsub("leader_|_bin", "", v))
  datbw$tmp <- datbw[[v]]
  if (sum(!is.na(datbw[[v]])) < 500) next
  mod <- feols(tmp ~ days_withdrawal * after_withdrawal | COUNTRYNEW[age] + COUNTRYNEW^educ + COUNTRYNEW^gender2, 
               cluster = ~ COUNTRYNEW,
               weights = ~HHWEIGHT2,
               data = datbw)
  
  res[[v]] <- tidy(mod) %>% 
    filter(term == "after_withdrawal") %>% 
    mutate(country = country) %>% 
    mutate(n = nobs(mod))
  
  datbw$tmp <- NULL
  
  cat(sprintf("Number of observations for %s: %s", country, nobs(mod)), "\n")
}
res <- bind_rows(res) %>% 
  mutate(country = case_when(
    country == "Usa" ~ "USA", 
    country == "Eu" ~ "EU",
    TRUE ~ country))


ggplot(res) + 
  aes(x = estimate, 
      xmin = estimate - 1.96 * std.error,
      xmax = estimate + 1.96 * std.error,
      y = fct_reorder(country, estimate),
      label = sprintf("N = %s", prettyNum(n, big.mark = ","))) + 
  geom_vline(xintercept = 0, lty = 3) + 
  geom_text(nudge_y = .25, hjust = 0) + 
  
  geom_point() + 
  geom_errorbarh(height = 0) + 
  labs(x = "Change in Approval of Country's Leader", y = "Country")
ggsave("Figures/RD-placebo-countries.pdf", width=6, height = 4)




## Russia and China tests of equality w/ US ------------------------------
zsum_base <- feols(
  c(leader_usa_bin, leader_china_bin, leader_russia_bin) ~ days_withdrawal * after_withdrawal | COUNTRYNEW,
  cluster = ~ COUNTRYNEW,
  weights = ~ HHWEIGHT2,
  data = datbw
)
zsum_full <- feols(
  c(leader_usa_bin, leader_china_bin, leader_russia_bin) ~ days_withdrawal * after_withdrawal |
    COUNTRYNEW[age] + COUNTRYNEW ^ educ + COUNTRYNEW ^ gender2,
  cluster = ~ COUNTRYNEW,
  weights = ~ HHWEIGHT2,
  data = datbw
)


# test USA = china
usa_china <- sur_hypotheses(zsum_full, 
                            hypothesis = "after_withdrawal_1 = after_withdrawal_2")
# test USA = russia
usa_russia <- sur_hypotheses(zsum_full,
                            hypothesis = "after_withdrawal_1 = after_withdrawal_3")


res %>% 
  filter(country %in% c("USA", "Russia", "China")) %>% 
  ggplot() + 
  aes(x = estimate, 
      xmin = estimate - 1.96 * std.error,
      xmax = estimate + 1.96 * std.error,
      y = fct_reorder(country, -1 *estimate),
      label = sprintf("N = %s", prettyNum(n, big.mark = ","))) + 
  geom_vline(xintercept = 0, lty = 3) + 
  geom_text(nudge_y = .25, hjust = .5) + 
  geom_point() + 
  geom_errorbarh(height = 0) + 
  geom_signif(comparisons = list(c("USA", "China"),
                                 c("USA", "Russia")),
              orientation = "y",
              y_position = c(.035, .025),
              tip_length = .01,
              annotations = c(paste0("p =", round(usa_china$p.value,3)),
                              paste0("p =", round(usa_russia$p.value,3)))) +
  
  scale_x_continuous(breaks = seq(-.1, .1, .01),
                     labels = ~ paste0(.x*100)) + 
  labs(x = "Percentage Point Change in\nApproval of Country's Leader", y = NULL)
ggsave("Figures/USA-china-russia.pdf", width=6, height = 4)






## Varying BW for Russia/China approval ---------------------------------
bws <- 10:100

### Russia ###
results_russia <- list()
i = 1
for (bw in bws){
  tmpdat <- dat %>% 
    filter(!is.na(leader_russia), 
           date_interview >= treatment_date - bw,
           date_interview <  treatment_date + bw) 
  mod1 <- feols(leader_russia_bin ~ days_withdrawal * after_withdrawal | COUNTRYNEW, 
                cluster = ~ COUNTRYNEW,
                weights = ~HHWEIGHT2,
                data = tmpdat)
  
  # use i() syntax for age covariate w/ varying slopes rather than COUNTRYNEW[age] 
  # due to convergence issues with the latter.
  mod2 <- feols(leader_russia_bin ~ days_withdrawal * after_withdrawal + i(COUNTRYNEW, age) |
                  COUNTRYNEW + COUNTRYNEW^educ + COUNTRYNEW^gender2, 
                cluster = ~ COUNTRYNEW,
                weights = ~HHWEIGHT2,
                data = tmpdat)
  
  
  out <- bind_rows(tidy(mod1) %>% mutate(model = "basic"),
                   tidy(mod2) %>% mutate(model = "controls")) %>% 
    filter(term == "after_withdrawal") %>% 
    mutate(bw = bw) 
  
  results_russia[[i]] <- out
  i = i+1
}
out_russia <- bind_rows(results_russia) %>% 
  mutate(target = "Russia")





### China ###
results_china <- list()
i = 1
for (bw in bws){
  tmpdat <- dat %>% 
    filter(!is.na(leader_china), 
           date_interview >= treatment_date - bw,
           date_interview <  treatment_date + bw) 
  mod1 <- feols(leader_china_bin ~ days_withdrawal * after_withdrawal | COUNTRYNEW, 
                cluster = ~ COUNTRYNEW,
                weights = ~HHWEIGHT2,
                data = tmpdat)
  
  # use i() syntax for age covariate w/ varying slopes rather than COUNTRYNEW[age] 
  # due to convergence issues with the latter.
  mod2 <- feols(leader_china_bin ~ days_withdrawal * after_withdrawal + i(COUNTRYNEW, age) |
                  COUNTRYNEW + COUNTRYNEW^educ + COUNTRYNEW^gender2, 
                cluster = ~ COUNTRYNEW,
                weights = ~HHWEIGHT2,
                data = tmpdat)
  
  
  out <- bind_rows(tidy(mod1) %>% mutate(model = "basic"),
                   tidy(mod2) %>% mutate(model = "controls")) %>% 
    filter(term == "after_withdrawal") %>% 
    mutate(bw = bw) 
  
  results_china[[i]] <- out
  i = i+1
}
out_china <- bind_rows(results_china) %>% 
  mutate(target = "China")




## Combine and plot
out_rc <- bind_rows(out_russia, out_china)
ggplot(out_rc) + 
  aes(x = bw, 
      y = estimate, 
      ymin = estimate - 1.96 * std.error, 
      ymax = estimate + 1.96 * std.error,
      colour = model, 
      fill = model) + 
  geom_hline(yintercept = 0, lty = 3) +
  geom_ribbon(alpha = .1) + 
  geom_point() + 
  scale_fill_manual(name = NULL, labels = c("basic", "w/ controls"), values = c("black", "blue")) + 
  scale_colour_manual(name = NULL, labels = c("basic", "w/ controls"), values = c("black", "blue")) + 
  labs(x = "Bandwidth (Days)", y = "RD Estimate") +
  theme(legend.position = c(.75, .15)) + 
  facet_wrap(~ target)
ggsave("Figures/RD-bandwith-plot-china-russia.pdf", width = 6, height = 4)





# Benchmarking the effect -------------------------------------------------

## Effect of Biden inauguration using same specifications (different data sadly)
bidendat <- alldat %>% 
  mutate(days_since_inaug = date_interview - as_date("2021-01-21"),
         biden = as.numeric(days_since_inaug >= 0)) %>% 
  filter(abs(days_since_inaug) <= 30) %>% 
  filter(!is.na(leader_usa))


## Get Sample Sizes by Country X Before/After
## Keep countries w/ at least 25 interviews both before/after
sampsizes <- bidendat %>% 
  count(COUNTRYNEW, factor(biden), .drop = FALSE) %>% 
  group_by(COUNTRYNEW) %>% 
  mutate(min = min(n)) %>% 
  filter(min >= 25)
tokeep <- unique(sampsizes$COUNTRYNEW)
bidendat <- bidendat %>% 
  filter(COUNTRYNEW %in% tokeep)
sampsizes <- sampsizes %>% 
  filter(COUNTRYNEW %in% bidendat$COUNTRYNEW)

# unfortunately only 4 countries in the field at the time of biden inaug
sampsizes

base <- feols(leader_usa_bin ~ biden,
              cluster = ~ COUNTRYNEW,
              weights = ~ HHWEIGHT2,
              data = bidendat)
basefe <- feols(leader_usa_bin ~ biden | COUNTRYNEW,
                cluster = ~ COUNTRYNEW,
                weights = ~ HHWEIGHT2,
                data = bidendat)
basefecov <- feols(leader_usa_bin ~ biden | COUNTRYNEW[age] + 
                     COUNTRYNEW^educ + COUNTRYNEW^gender2,
                   cluster = ~ COUNTRYNEW,
                   weights = ~ HHWEIGHT2,
                   data = bidendat)
rd <- feols(leader_usa_bin ~ biden * days_since_inaug, 
            cluster = ~ COUNTRYNEW,
            weights = ~HHWEIGHT2,
            data = bidendat)
rdfe <- feols(leader_usa_bin ~ biden * days_since_inaug | COUNTRYNEW, 
              cluster = ~ COUNTRYNEW,
              weights = ~HHWEIGHT2,
              data = bidendat)
rdfecov <- feols(leader_usa_bin ~ biden * days_since_inaug | COUNTRYNEW[age] + 
                   COUNTRYNEW^educ + COUNTRYNEW^gender2, 
                 cluster = ~ COUNTRYNEW,
                 weights = ~HHWEIGHT2,
                 data = bidendat)



modelsummary(list(base, basefe, basefecov, rd, rdfe, rdfecov), 
             stars = c("$^*$" = .1, "$^{**}$" = .05, "$^{***}$" = .01), 
             coef_map = c("biden" = "Biden Inaugurated",
                          "days_since_inaug" = "Time",
                          "biden:days_since_inaug" = "Time $\\times$ Biden Inaugurated"),
             booktabs = TRUE, 
             gof_omit = c("Std.Errors|AIC|BIC|R2 Adj."),
             output = "latex_tabular",
             escape = FALSE,
             gof_map = tribble(~ raw, ~ clean, ~ fmt, 
                               "nobs", "$N$", function(.x) prettyNum(.x, big.mark = ","),
                               "FE: COUNTRYNEW", "Country FE", modelsummary::fmt_statistic),
             add_rows = tribble(~ term, ~ base, ~ basefe, ~ basefecov, ~ rd, ~ rdfe, ~ rdfecov,
                                "Ind. Covariates $\\times$ Country FE", "", "", "X", "", "", "X")
) %>% 
  format_tt(replace = list("$\\checkmark$" = "X")) %>% 
  save_tt("Tables/biden-inaug.tex", overwrite = TRUE)






# Balance tests -----------------------------------------------------------

# Balance test: predict whether respondent was interviewed before/after as
# function of baseline covariates

balance <- feols(
  after_withdrawal ~ age + male + educ + hhsize + inc_quint_2 +
    inc_quint_3 + inc_quint_4 + inc_quint_5,
  cluster = ~ COUNTRYNEW,
  weights = ~ HHWEIGHT2,
  data = datbw
)

balance_fe <- feols(
  after_withdrawal ~ age + male + educ + hhsize + inc_quint_2 +
    inc_quint_3 + inc_quint_4 + inc_quint_5  |
    COUNTRYNEW,
  cluster = ~ COUNTRYNEW,
  weights = ~ HHWEIGHT2,
  data = datbw
)






modelsummary(list(balance, balance_fe), 
             stars = c("$^*$" = .1, "$^{**}$" = .05, "$^{***}$" = .01), 
             coef_map = c("age" = "Age",
                          "male" = "Male",
                          "educSecondary school" = "Secondary educ.",
                          "educPostsecondary school" = "Post-secondary educ.",
                          "hhsize" = "Household size",
                          "inc_quint_2" = "Income quintile: 2nd",
                          "inc_quint_3" = "Income quintile: 3rd",
                          "inc_quint_4" = "Income quintile: 4th",
                          "inc_quint_5" = "Income quintile: 5th"),
             booktabs = TRUE, 
             gof_omit = c("Std.Errors|AIC|BIC|R2 Adj."),
             output = "latex_tabular",
             escape = FALSE,
             gof_map = tribble(~ raw, ~ clean, ~ fmt, 
                               "nobs", "$N$", function(.x) prettyNum(.x, big.mark = ","),
                               "FE: COUNTRYNEW", "Country FE", modelsummary::fmt_statistic)) %>% 
  format_tt(replace = list("$\\checkmark$" = "X")) %>% 
  save_tt("Tables/balance-regressions.tex", overwrite = TRUE)







# rdrobust ----------------------------------------------------------------

message("Running rdrobust models...")

rddat <- dat %>% 
  filter(!is.na(leader_usa), abs(days_withdrawal) <= 90)
country.fe <- model.matrix( ~ rddat$COUNTRYNEW + 0)

rddat_covs <- rddat %>% 
  drop_na(age, educ, gender2)
full.covs <-  model.matrix( ~ COUNTRYNEW * (educ + age + gender2) +  0,  data = rddat_covs)

rdrob <- rdrobust(y = rddat$leader_usa_bin, x = rddat$days_withdrawal, 
                  covs = country.fe, weights = rddat$HHWEIGHT2)

rdrob_cov <- rdrobust(y = rddat_covs$leader_usa_bin, x = rddat_covs$days_withdrawal, 
                      covs = full.covs, weights = rddat_covs$HHWEIGHT2)

# function to generate table from RD robust
rdrobust_table <- function(model, file = NULL) {
  coef <- model$coef
  se   <- model$se
  ci   <- apply(model$ci, 1, \(x) sprintf("[%s]", paste0(round(x, 3), collapse = ", ")))
  bw   <- model$bws[1,1]
  pval <- model$pv
  tab <- tibble(
    Type = rownames(coef),
    Estimate = round(coef[, 1], 3),
    `Std. Error` = round(se[, 1], 4),
    `95\\% CI` = ci,
    `p-value` = round(pval[, 1], 3),
    Bandwidth = bw
  )
  # turn into latex tabular environment
  out <- tab %>% 
    xtable::xtable(digits = 3) %>% 
    print(include.rownames = FALSE, 
          sanitize.colnames.function = function(x) {paste('{\\textbf{',x,'}}', sep = '')},
          booktabs = TRUE, 
          floating = FALSE,
          escape = FALSE,
          add.to.row = list(pos = list(-1), 
                            command = c("\\toprule")),
          hline.after = c(0, nrow(tab))) 
  if (!is.null(file)) {
    writeLines(out, con = file)
  }
}

rdrobust_table(rdrob, file = "Tables/rdrobust-fe.tex")
rdrobust_table(rdrob_cov, file = "Tables/rdrobust-covs.tex")

cat(round(rdrob$bws[1,1], 1), file = "Tables/cct-bw.txt")
cat(round(rdrob_cov$bws[1,1], 1), file = "Tables/cct-bw-covs.txt")


# save models
saveRDS(list(rdrob = rdrob, rdrob_cov = rdrob_cov), file = "Data/rdrobust-fits.rds")
